library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.0 ✔ purrr 0.3.2
## ✔ tibble 2.1.1 ✔ dplyr 0.8.3
## ✔ tidyr 1.0.0 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(glue)
##
## Attaching package: 'glue'
## The following object is masked from 'package:dplyr':
##
## collapse
library(tidybayes)
library(cowplot)
##
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
## default ggplot2 theme anymore. To recover the previous
## behavior, execute:
## theme_set(theme_cowplot())
## ********************************************************
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
library(GGally)
##
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
##
## nasa
library(knitr)
theme_set(theme_cowplot())
load_existing_sim <- TRUE
Continuing from m1+m2, the most noticeable problem with the initial priors in m1 was that too many extremes values were being observed at the subject level, both in terms of cricSD and pMem. I think m2 achieved good priors on pMem, but I think in addition to narrowing variance priors on circSD I need to shift the mean prior.
$$ \begin{aligned}
mu__0 Normal(4, 0.5) ~~ …&to… ~~ mu__0 Normal(3.8, 0.5) \ \ sigma_0 Normal^+(0, 0.5) ~~ …&to… ~~ sigma_0 Normal^+(0, 0.25) \ sigma_0 Normal^+(0, 1.5) ~~ …&to… ~~ sigma_0 Normal^+(0, 0.5) \ \ sigma_Normal^+(0, 0.5) ~~ …&to… ~~ sigma_Normal^+(0, 0.25) \ sigma_Normal^+(0, 1) ~~ …&to… ~~ sigma_Normal^+(0, 0.5) \end{aligned} $$
Written down the model is:
\[ \begin{aligned} \mathrm{-likelihood-} \\ error_i &\sim (pMem_i)*VM(0, \kappa_i) + (1 - pMem_i)*Unif(-\pi,\pi) \\ \\ \mathrm{-param~transformation-} \\ \kappa_i &= sd2k(circ\_sd_i) \\ \\ \mathrm{-linear~model-} \\ circ\_sd_i &= exp(\alpha_{0,SUBJ[i]} + \alpha_{\Delta,SUBJ[i]} * postCond) ~~ \mathrm{-log~link~on~circSD-} \\ pMem_i &= inv\_logit(\beta_{0,SUBJ[i]} + \beta_{\Delta,SUBJ[i]} * postCond) ~~ \mathrm{-logit~link~on~pMem-} \\ \\ \mathrm{-priors:~all~independent-}\\ \alpha_{0,SUBJ[...]} &\sim Normal(mu\_\alpha_0, sigma\_\alpha_0) \\ \alpha_{\Delta,SUBJ[...]} &\sim Normal(mu\_\alpha_{\Delta}, sigma\_\alpha_{\Delta}) \\ \beta_{0,SUBJ[...]} &\sim Normal(mu\_\beta_0, sigma\_\beta_0) \\ \beta_{\Delta,SUBJ[...]} &\sim Normal(mu\_\beta_{\Delta}, sigma\_\beta_{\Delta}) \\ \\ mu\_\alpha_0 &\sim Normal(3.8, 0.5)\\ sigma\_\alpha_0 &\sim Normal^+(0, 0.25) \\ mu\_\alpha_{\Delta} &\sim Normal(0, 0.5)\\ sigma\_\alpha_{\Delta} &\sim Normal^+(0, 0.25)\\ mu\_\beta_0 &\sim Normal(0, 1.5)\\ sigma\_\beta_0 &\sim Normal^+(0, 0.5)\\ mu\_\beta_{\Delta} &\sim Normal(0, 1)\\ sigma\_\beta_{\Delta} &\sim Normal^+(0, 0.5)\\ \end{aligned} \]
I am considering the color stimulation data. samples sizes:
obs_data <- read_csv('../data/stimulation_obvs.csv')
## Parsed with column specification:
## cols(
## subj = col_double(),
## subj_index = col_double(),
## stimulation = col_double(),
## error = col_double()
## )
#summarize subj num, obs count, and obs condition split
(subj_summary <- obs_data %>%
group_by(subj_index) %>%
summarise(n_obs = n(), frac_stim = mean(stimulation)))
## # A tibble: 2 x 3
## subj_index n_obs frac_stim
## <dbl> <int> <dbl>
## 1 1 252 0.5
## 2 2 252 0.5
2 subjs with 252 observations each, 126 per condition.
9/12 - I am thinking I should simulate varying effects using the actual samples sizes I have. I also think that this shouldnt matter (at least for a model like this). I’ll also try with a larger number of subjects.
nsubj_sim <- 15
nobs_per_cond_sim <- 126
functions
inv_logit_scalar <- function(x){
return(exp(x)/(exp(x) + 1))
}
inv_logit <- Vectorize(inv_logit_scalar)
sd2k <- function(sd_input) {
#sd_input needs to be radians
#original matlab code
# R <- exp(-S.^2/2);
# K = 1./(R.^3 - 4 * R.^2 + 3 * R);
# K(R < 0.85) = -0.4 + 1.39 * R(R < 0.85) + 0.43./(1 - R(R < 0.85));
# K(R < 0.53) = 2 * R(R < 0.53) + R(R < 0.53).^3 + (5 * R(R < 0.53).^5)/6;
R <- exp(-(sd_input^2)/2)
K <- 1/((R^3) - (4 * R^2) + (3*R))
if(R < 0.85){
K <- -0.4 + (1.39 * R) + (0.43/(1-R))
}
if(R < 0.53){
K <- (2 * R) + R^3 + (5 * R^5)/6
}
return(K)
}
sd2k_vec <- Vectorize(sd2k)
# this function takes a single set of parameters and covariates and simlates an observation from the mixture model
simulateData_likelihood <- function(pMem,
k,
nobs = 1){
# which part of the mixture does this obs come from
memFlip <- rbernoulli(nobs, pMem)
obs <- memFlip * ( Rfast::rvonmises(nobs, pi, k) - pi) + (1 - memFlip) * runif(nobs, -pi, pi)
obs <- as.numeric(obs)
obs_tibble <- tibble(obs_radian = obs, obs_degree = pracma::rad2deg(obs))
return(obs_tibble)
}
simulateData_subj <- function(subj_alpha0, subj_alphaD,
subj_beta0, subj_betaD,
nobs_per_condition){
conditions <- c(0,1)
stimulation <- rep(conditions, each = nobs_per_condition)
params <- list(pMem = inv_logit(subj_beta0 + (subj_betaD * stimulation)),
k = sd2k_vec(
pracma::deg2rad(
exp(subj_alpha0 + (subj_alphaD * stimulation)))))
subj_obs <- pmap_dfr(params, simulateData_likelihood) %>%
cbind(stimulation)
return(subj_obs)
}
# this function will take a single set of simulated parameter draws and then simulate an implied data set
simulateData <- function(sim_num,
alpha0_mu, alpha0_sigma,
alphaD_mu, alphaD_sigma,
beta0_mu, beta0_sigma,
betaD_mu, betaD_sigma,
nsubj = 100, nobs_per_condition = 100){
print(glue("simulation: {sim_num}"))
# generate lower level subject-specific parameters
sim_subj_alpha0 <- rnorm(nsubj, alpha0_mu, alpha0_sigma)
sim_subj_alphaD <- rnorm(nsubj, alphaD_mu, alphaD_sigma)
sim_subj_beta0 <- rnorm(nsubj, beta0_mu, beta0_sigma)
sim_subj_betaD <- rnorm(nsubj, betaD_mu, betaD_sigma)
sim_data <- tibble(
subj = 1:nsubj,
subj_alpha0 = sim_subj_alpha0,
subj_alphaD = sim_subj_alphaD,
subj_beta0 = sim_subj_beta0,
subj_betaD = sim_subj_betaD,
nobs_per_condition = nobs_per_condition
)
sim_data <- sim_data %>%
mutate(subj_obs = pmap(sim_data %>% select(-subj), simulateData_subj))
return(sim_data)
}
run simulate
found_existing_sim <- FALSE
if (file.exists("sim_datasets.rds")){
found_existing_sim <- TRUE
}
if (load_existing_sim == FALSE || found_existing_sim == FALSE){
nsim_datasets <- 1200
sim_priors <- tibble(
sim_num = 1:nsim_datasets,
alpha0_mu = rnorm(nsim_datasets, 3.8, 0.5),
alpha0_sigma = abs(rnorm(nsim_datasets, 0, 0.25)),
alphaD_mu = rnorm(nsim_datasets, 0, 0.5),
alphaD_sigma = abs(rnorm(nsim_datasets, 0, 0.25)),
beta0_mu = rnorm(nsim_datasets, 0, 1.5),
beta0_sigma = abs(rnorm(nsim_datasets, 0, 0.5)),
betaD_mu = rnorm(nsim_datasets, 0, 1),
betaD_sigma = abs(rnorm(nsim_datasets, 0, 0.5)),
nsubj = nsubj_sim,
nobs_per_cond = nobs_per_cond_sim
)
sim_datasets <- sim_priors %>%
mutate(dataset = pmap(sim_priors, simulateData))
saveRDS(sim_datasets, file = "sim_datasets.rds")
}else{
sim_datasets <- readRDS("sim_datasets.rds")
}
## simulation: 1
## simulation: 2
## simulation: 3
## simulation: 4
## simulation: 5
## simulation: 6
## simulation: 7
## simulation: 8
## simulation: 9
## simulation: 10
## simulation: 11
## simulation: 12
## simulation: 13
## simulation: 14
## simulation: 15
## simulation: 16
## simulation: 17
## simulation: 18
## simulation: 19
## simulation: 20
## simulation: 21
## simulation: 22
## simulation: 23
## simulation: 24
## simulation: 25
## simulation: 26
## simulation: 27
## simulation: 28
## simulation: 29
## simulation: 30
## simulation: 31
## simulation: 32
## simulation: 33
## simulation: 34
## simulation: 35
## simulation: 36
## simulation: 37
## simulation: 38
## simulation: 39
## simulation: 40
## simulation: 41
## simulation: 42
## simulation: 43
## simulation: 44
## simulation: 45
## simulation: 46
## simulation: 47
## simulation: 48
## simulation: 49
## simulation: 50
## simulation: 51
## simulation: 52
## simulation: 53
## simulation: 54
## simulation: 55
## simulation: 56
## simulation: 57
## simulation: 58
## simulation: 59
## simulation: 60
## simulation: 61
## simulation: 62
## simulation: 63
## simulation: 64
## simulation: 65
## simulation: 66
## simulation: 67
## simulation: 68
## simulation: 69
## simulation: 70
## simulation: 71
## simulation: 72
## simulation: 73
## simulation: 74
## simulation: 75
## simulation: 76
## simulation: 77
## simulation: 78
## simulation: 79
## simulation: 80
## simulation: 81
## simulation: 82
## simulation: 83
## simulation: 84
## simulation: 85
## simulation: 86
## simulation: 87
## simulation: 88
## simulation: 89
## simulation: 90
## simulation: 91
## simulation: 92
## simulation: 93
## simulation: 94
## simulation: 95
## simulation: 96
## simulation: 97
## simulation: 98
## simulation: 99
## simulation: 100
## simulation: 101
## simulation: 102
## simulation: 103
## simulation: 104
## simulation: 105
## simulation: 106
## simulation: 107
## simulation: 108
## simulation: 109
## simulation: 110
## simulation: 111
## simulation: 112
## simulation: 113
## simulation: 114
## simulation: 115
## simulation: 116
## simulation: 117
## simulation: 118
## simulation: 119
## simulation: 120
## simulation: 121
## simulation: 122
## simulation: 123
## simulation: 124
## simulation: 125
## simulation: 126
## simulation: 127
## simulation: 128
## simulation: 129
## simulation: 130
## simulation: 131
## simulation: 132
## simulation: 133
## simulation: 134
## simulation: 135
## simulation: 136
## simulation: 137
## simulation: 138
## simulation: 139
## simulation: 140
## simulation: 141
## simulation: 142
## simulation: 143
## simulation: 144
## simulation: 145
## simulation: 146
## simulation: 147
## simulation: 148
## simulation: 149
## simulation: 150
## simulation: 151
## simulation: 152
## simulation: 153
## simulation: 154
## simulation: 155
## simulation: 156
## simulation: 157
## simulation: 158
## simulation: 159
## simulation: 160
## simulation: 161
## simulation: 162
## simulation: 163
## simulation: 164
## simulation: 165
## simulation: 166
## simulation: 167
## simulation: 168
## simulation: 169
## simulation: 170
## simulation: 171
## simulation: 172
## simulation: 173
## simulation: 174
## simulation: 175
## simulation: 176
## simulation: 177
## simulation: 178
## simulation: 179
## simulation: 180
## simulation: 181
## simulation: 182
## simulation: 183
## simulation: 184
## simulation: 185
## simulation: 186
## simulation: 187
## simulation: 188
## simulation: 189
## simulation: 190
## simulation: 191
## simulation: 192
## simulation: 193
## simulation: 194
## simulation: 195
## simulation: 196
## simulation: 197
## simulation: 198
## simulation: 199
## simulation: 200
## simulation: 201
## simulation: 202
## simulation: 203
## simulation: 204
## simulation: 205
## simulation: 206
## simulation: 207
## simulation: 208
## simulation: 209
## simulation: 210
## simulation: 211
## simulation: 212
## simulation: 213
## simulation: 214
## simulation: 215
## simulation: 216
## simulation: 217
## simulation: 218
## simulation: 219
## simulation: 220
## simulation: 221
## simulation: 222
## simulation: 223
## simulation: 224
## simulation: 225
## simulation: 226
## simulation: 227
## simulation: 228
## simulation: 229
## simulation: 230
## simulation: 231
## simulation: 232
## simulation: 233
## simulation: 234
## simulation: 235
## simulation: 236
## simulation: 237
## simulation: 238
## simulation: 239
## simulation: 240
## simulation: 241
## simulation: 242
## simulation: 243
## simulation: 244
## simulation: 245
## simulation: 246
## simulation: 247
## simulation: 248
## simulation: 249
## simulation: 250
## simulation: 251
## simulation: 252
## simulation: 253
## simulation: 254
## simulation: 255
## simulation: 256
## simulation: 257
## simulation: 258
## simulation: 259
## simulation: 260
## simulation: 261
## simulation: 262
## simulation: 263
## simulation: 264
## simulation: 265
## simulation: 266
## simulation: 267
## simulation: 268
## simulation: 269
## simulation: 270
## simulation: 271
## simulation: 272
## simulation: 273
## simulation: 274
## simulation: 275
## simulation: 276
## simulation: 277
## simulation: 278
## simulation: 279
## simulation: 280
## simulation: 281
## simulation: 282
## simulation: 283
## simulation: 284
## simulation: 285
## simulation: 286
## simulation: 287
## simulation: 288
## simulation: 289
## simulation: 290
## simulation: 291
## simulation: 292
## simulation: 293
## simulation: 294
## simulation: 295
## simulation: 296
## simulation: 297
## simulation: 298
## simulation: 299
## simulation: 300
## simulation: 301
## simulation: 302
## simulation: 303
## simulation: 304
## simulation: 305
## simulation: 306
## simulation: 307
## simulation: 308
## simulation: 309
## simulation: 310
## simulation: 311
## simulation: 312
## simulation: 313
## simulation: 314
## simulation: 315
## simulation: 316
## simulation: 317
## simulation: 318
## simulation: 319
## simulation: 320
## simulation: 321
## simulation: 322
## simulation: 323
## simulation: 324
## simulation: 325
## simulation: 326
## simulation: 327
## simulation: 328
## simulation: 329
## simulation: 330
## simulation: 331
## simulation: 332
## simulation: 333
## simulation: 334
## simulation: 335
## simulation: 336
## simulation: 337
## simulation: 338
## simulation: 339
## simulation: 340
## simulation: 341
## simulation: 342
## simulation: 343
## simulation: 344
## simulation: 345
## simulation: 346
## simulation: 347
## simulation: 348
## simulation: 349
## simulation: 350
## simulation: 351
## simulation: 352
## simulation: 353
## simulation: 354
## simulation: 355
## simulation: 356
## simulation: 357
## simulation: 358
## simulation: 359
## simulation: 360
## simulation: 361
## simulation: 362
## simulation: 363
## simulation: 364
## simulation: 365
## simulation: 366
## simulation: 367
## simulation: 368
## simulation: 369
## simulation: 370
## simulation: 371
## simulation: 372
## simulation: 373
## simulation: 374
## simulation: 375
## simulation: 376
## simulation: 377
## simulation: 378
## simulation: 379
## simulation: 380
## simulation: 381
## simulation: 382
## simulation: 383
## simulation: 384
## simulation: 385
## simulation: 386
## simulation: 387
## simulation: 388
## simulation: 389
## simulation: 390
## simulation: 391
## simulation: 392
## simulation: 393
## simulation: 394
## simulation: 395
## simulation: 396
## simulation: 397
## simulation: 398
## simulation: 399
## simulation: 400
## simulation: 401
## simulation: 402
## simulation: 403
## simulation: 404
## simulation: 405
## simulation: 406
## simulation: 407
## simulation: 408
## simulation: 409
## simulation: 410
## simulation: 411
## simulation: 412
## simulation: 413
## simulation: 414
## simulation: 415
## simulation: 416
## simulation: 417
## simulation: 418
## simulation: 419
## simulation: 420
## simulation: 421
## simulation: 422
## simulation: 423
## simulation: 424
## simulation: 425
## simulation: 426
## simulation: 427
## simulation: 428
## simulation: 429
## simulation: 430
## simulation: 431
## simulation: 432
## simulation: 433
## simulation: 434
## simulation: 435
## simulation: 436
## simulation: 437
## simulation: 438
## simulation: 439
## simulation: 440
## simulation: 441
## simulation: 442
## simulation: 443
## simulation: 444
## simulation: 445
## simulation: 446
## simulation: 447
## simulation: 448
## simulation: 449
## simulation: 450
## simulation: 451
## simulation: 452
## simulation: 453
## simulation: 454
## simulation: 455
## simulation: 456
## simulation: 457
## simulation: 458
## simulation: 459
## simulation: 460
## simulation: 461
## simulation: 462
## simulation: 463
## simulation: 464
## simulation: 465
## simulation: 466
## simulation: 467
## simulation: 468
## simulation: 469
## simulation: 470
## simulation: 471
## simulation: 472
## simulation: 473
## simulation: 474
## simulation: 475
## simulation: 476
## simulation: 477
## simulation: 478
## simulation: 479
## simulation: 480
## simulation: 481
## simulation: 482
## simulation: 483
## simulation: 484
## simulation: 485
## simulation: 486
## simulation: 487
## simulation: 488
## simulation: 489
## simulation: 490
## simulation: 491
## simulation: 492
## simulation: 493
## simulation: 494
## simulation: 495
## simulation: 496
## simulation: 497
## simulation: 498
## simulation: 499
## simulation: 500
## simulation: 501
## simulation: 502
## simulation: 503
## simulation: 504
## simulation: 505
## simulation: 506
## simulation: 507
## simulation: 508
## simulation: 509
## simulation: 510
## simulation: 511
## simulation: 512
## simulation: 513
## simulation: 514
## simulation: 515
## simulation: 516
## simulation: 517
## simulation: 518
## simulation: 519
## simulation: 520
## simulation: 521
## simulation: 522
## simulation: 523
## simulation: 524
## simulation: 525
## simulation: 526
## simulation: 527
## simulation: 528
## simulation: 529
## simulation: 530
## simulation: 531
## simulation: 532
## simulation: 533
## simulation: 534
## simulation: 535
## simulation: 536
## simulation: 537
## simulation: 538
## simulation: 539
## simulation: 540
## simulation: 541
## simulation: 542
## simulation: 543
## simulation: 544
## simulation: 545
## simulation: 546
## simulation: 547
## simulation: 548
## simulation: 549
## simulation: 550
## simulation: 551
## simulation: 552
## simulation: 553
## simulation: 554
## simulation: 555
## simulation: 556
## simulation: 557
## simulation: 558
## simulation: 559
## simulation: 560
## simulation: 561
## simulation: 562
## simulation: 563
## simulation: 564
## simulation: 565
## simulation: 566
## simulation: 567
## simulation: 568
## simulation: 569
## simulation: 570
## simulation: 571
## simulation: 572
## simulation: 573
## simulation: 574
## simulation: 575
## simulation: 576
## simulation: 577
## simulation: 578
## simulation: 579
## simulation: 580
## simulation: 581
## simulation: 582
## simulation: 583
## simulation: 584
## simulation: 585
## simulation: 586
## simulation: 587
## simulation: 588
## simulation: 589
## simulation: 590
## simulation: 591
## simulation: 592
## simulation: 593
## simulation: 594
## simulation: 595
## simulation: 596
## simulation: 597
## simulation: 598
## simulation: 599
## simulation: 600
## simulation: 601
## simulation: 602
## simulation: 603
## simulation: 604
## simulation: 605
## simulation: 606
## simulation: 607
## simulation: 608
## simulation: 609
## simulation: 610
## simulation: 611
## simulation: 612
## simulation: 613
## simulation: 614
## simulation: 615
## simulation: 616
## simulation: 617
## simulation: 618
## simulation: 619
## simulation: 620
## simulation: 621
## simulation: 622
## simulation: 623
## simulation: 624
## simulation: 625
## simulation: 626
## simulation: 627
## simulation: 628
## simulation: 629
## simulation: 630
## simulation: 631
## simulation: 632
## simulation: 633
## simulation: 634
## simulation: 635
## simulation: 636
## simulation: 637
## simulation: 638
## simulation: 639
## simulation: 640
## simulation: 641
## simulation: 642
## simulation: 643
## simulation: 644
## simulation: 645
## simulation: 646
## simulation: 647
## simulation: 648
## simulation: 649
## simulation: 650
## simulation: 651
## simulation: 652
## simulation: 653
## simulation: 654
## simulation: 655
## simulation: 656
## simulation: 657
## simulation: 658
## simulation: 659
## simulation: 660
## simulation: 661
## simulation: 662
## simulation: 663
## simulation: 664
## simulation: 665
## simulation: 666
## simulation: 667
## simulation: 668
## simulation: 669
## simulation: 670
## simulation: 671
## simulation: 672
## simulation: 673
## simulation: 674
## simulation: 675
## simulation: 676
## simulation: 677
## simulation: 678
## simulation: 679
## simulation: 680
## simulation: 681
## simulation: 682
## simulation: 683
## simulation: 684
## simulation: 685
## simulation: 686
## simulation: 687
## simulation: 688
## simulation: 689
## simulation: 690
## simulation: 691
## simulation: 692
## simulation: 693
## simulation: 694
## simulation: 695
## simulation: 696
## simulation: 697
## simulation: 698
## simulation: 699
## simulation: 700
## simulation: 701
## simulation: 702
## simulation: 703
## simulation: 704
## simulation: 705
## simulation: 706
## simulation: 707
## simulation: 708
## simulation: 709
## simulation: 710
## simulation: 711
## simulation: 712
## simulation: 713
## simulation: 714
## simulation: 715
## simulation: 716
## simulation: 717
## simulation: 718
## simulation: 719
## simulation: 720
## simulation: 721
## simulation: 722
## simulation: 723
## simulation: 724
## simulation: 725
## simulation: 726
## simulation: 727
## simulation: 728
## simulation: 729
## simulation: 730
## simulation: 731
## simulation: 732
## simulation: 733
## simulation: 734
## simulation: 735
## simulation: 736
## simulation: 737
## simulation: 738
## simulation: 739
## simulation: 740
## simulation: 741
## simulation: 742
## simulation: 743
## simulation: 744
## simulation: 745
## simulation: 746
## simulation: 747
## simulation: 748
## simulation: 749
## simulation: 750
## simulation: 751
## simulation: 752
## simulation: 753
## simulation: 754
## simulation: 755
## simulation: 756
## simulation: 757
## simulation: 758
## simulation: 759
## simulation: 760
## simulation: 761
## simulation: 762
## simulation: 763
## simulation: 764
## simulation: 765
## simulation: 766
## simulation: 767
## simulation: 768
## simulation: 769
## simulation: 770
## simulation: 771
## simulation: 772
## simulation: 773
## simulation: 774
## simulation: 775
## simulation: 776
## simulation: 777
## simulation: 778
## simulation: 779
## simulation: 780
## simulation: 781
## simulation: 782
## simulation: 783
## simulation: 784
## simulation: 785
## simulation: 786
## simulation: 787
## simulation: 788
## simulation: 789
## simulation: 790
## simulation: 791
## simulation: 792
## simulation: 793
## simulation: 794
## simulation: 795
## simulation: 796
## simulation: 797
## simulation: 798
## simulation: 799
## simulation: 800
## simulation: 801
## simulation: 802
## simulation: 803
## simulation: 804
## simulation: 805
## simulation: 806
## simulation: 807
## simulation: 808
## simulation: 809
## simulation: 810
## simulation: 811
## simulation: 812
## simulation: 813
## simulation: 814
## simulation: 815
## simulation: 816
## simulation: 817
## simulation: 818
## simulation: 819
## simulation: 820
## simulation: 821
## simulation: 822
## simulation: 823
## simulation: 824
## simulation: 825
## simulation: 826
## simulation: 827
## simulation: 828
## simulation: 829
## simulation: 830
## simulation: 831
## simulation: 832
## simulation: 833
## simulation: 834
## simulation: 835
## simulation: 836
## simulation: 837
## simulation: 838
## simulation: 839
## simulation: 840
## simulation: 841
## simulation: 842
## simulation: 843
## simulation: 844
## simulation: 845
## simulation: 846
## simulation: 847
## simulation: 848
## simulation: 849
## simulation: 850
## simulation: 851
## simulation: 852
## simulation: 853
## simulation: 854
## simulation: 855
## simulation: 856
## simulation: 857
## simulation: 858
## simulation: 859
## simulation: 860
## simulation: 861
## simulation: 862
## simulation: 863
## simulation: 864
## simulation: 865
## simulation: 866
## simulation: 867
## simulation: 868
## simulation: 869
## simulation: 870
## simulation: 871
## simulation: 872
## simulation: 873
## simulation: 874
## simulation: 875
## simulation: 876
## simulation: 877
## simulation: 878
## simulation: 879
## simulation: 880
## simulation: 881
## simulation: 882
## simulation: 883
## simulation: 884
## simulation: 885
## simulation: 886
## simulation: 887
## simulation: 888
## simulation: 889
## simulation: 890
## simulation: 891
## simulation: 892
## simulation: 893
## simulation: 894
## simulation: 895
## simulation: 896
## simulation: 897
## simulation: 898
## simulation: 899
## simulation: 900
## simulation: 901
## simulation: 902
## simulation: 903
## simulation: 904
## simulation: 905
## simulation: 906
## simulation: 907
## simulation: 908
## simulation: 909
## simulation: 910
## simulation: 911
## simulation: 912
## simulation: 913
## simulation: 914
## simulation: 915
## simulation: 916
## simulation: 917
## simulation: 918
## simulation: 919
## simulation: 920
## simulation: 921
## simulation: 922
## simulation: 923
## simulation: 924
## simulation: 925
## simulation: 926
## simulation: 927
## simulation: 928
## simulation: 929
## simulation: 930
## simulation: 931
## simulation: 932
## simulation: 933
## simulation: 934
## simulation: 935
## simulation: 936
## simulation: 937
## simulation: 938
## simulation: 939
## simulation: 940
## simulation: 941
## simulation: 942
## simulation: 943
## simulation: 944
## simulation: 945
## simulation: 946
## simulation: 947
## simulation: 948
## simulation: 949
## simulation: 950
## simulation: 951
## simulation: 952
## simulation: 953
## simulation: 954
## simulation: 955
## simulation: 956
## simulation: 957
## simulation: 958
## simulation: 959
## simulation: 960
## simulation: 961
## simulation: 962
## simulation: 963
## simulation: 964
## simulation: 965
## simulation: 966
## simulation: 967
## simulation: 968
## simulation: 969
## simulation: 970
## simulation: 971
## simulation: 972
## simulation: 973
## simulation: 974
## simulation: 975
## simulation: 976
## simulation: 977
## simulation: 978
## simulation: 979
## simulation: 980
## simulation: 981
## simulation: 982
## simulation: 983
## simulation: 984
## simulation: 985
## simulation: 986
## simulation: 987
## simulation: 988
## simulation: 989
## simulation: 990
## simulation: 991
## simulation: 992
## simulation: 993
## simulation: 994
## simulation: 995
## simulation: 996
## simulation: 997
## simulation: 998
## simulation: 999
## simulation: 1000
## simulation: 1001
## simulation: 1002
## simulation: 1003
## simulation: 1004
## simulation: 1005
## simulation: 1006
## simulation: 1007
## simulation: 1008
## simulation: 1009
## simulation: 1010
## simulation: 1011
## simulation: 1012
## simulation: 1013
## simulation: 1014
## simulation: 1015
## simulation: 1016
## simulation: 1017
## simulation: 1018
## simulation: 1019
## simulation: 1020
## simulation: 1021
## simulation: 1022
## simulation: 1023
## simulation: 1024
## simulation: 1025
## simulation: 1026
## simulation: 1027
## simulation: 1028
## simulation: 1029
## simulation: 1030
## simulation: 1031
## simulation: 1032
## simulation: 1033
## simulation: 1034
## simulation: 1035
## simulation: 1036
## simulation: 1037
## simulation: 1038
## simulation: 1039
## simulation: 1040
## simulation: 1041
## simulation: 1042
## simulation: 1043
## simulation: 1044
## simulation: 1045
## simulation: 1046
## simulation: 1047
## simulation: 1048
## simulation: 1049
## simulation: 1050
## simulation: 1051
## simulation: 1052
## simulation: 1053
## simulation: 1054
## simulation: 1055
## simulation: 1056
## simulation: 1057
## simulation: 1058
## simulation: 1059
## simulation: 1060
## simulation: 1061
## simulation: 1062
## simulation: 1063
## simulation: 1064
## simulation: 1065
## simulation: 1066
## simulation: 1067
## simulation: 1068
## simulation: 1069
## simulation: 1070
## simulation: 1071
## simulation: 1072
## simulation: 1073
## simulation: 1074
## simulation: 1075
## simulation: 1076
## simulation: 1077
## simulation: 1078
## simulation: 1079
## simulation: 1080
## simulation: 1081
## simulation: 1082
## simulation: 1083
## simulation: 1084
## simulation: 1085
## simulation: 1086
## simulation: 1087
## simulation: 1088
## simulation: 1089
## simulation: 1090
## simulation: 1091
## simulation: 1092
## simulation: 1093
## simulation: 1094
## simulation: 1095
## simulation: 1096
## simulation: 1097
## simulation: 1098
## simulation: 1099
## simulation: 1100
## simulation: 1101
## simulation: 1102
## simulation: 1103
## simulation: 1104
## simulation: 1105
## simulation: 1106
## simulation: 1107
## simulation: 1108
## simulation: 1109
## simulation: 1110
## simulation: 1111
## simulation: 1112
## simulation: 1113
## simulation: 1114
## simulation: 1115
## simulation: 1116
## simulation: 1117
## simulation: 1118
## simulation: 1119
## simulation: 1120
## simulation: 1121
## simulation: 1122
## simulation: 1123
## simulation: 1124
## simulation: 1125
## simulation: 1126
## simulation: 1127
## simulation: 1128
## simulation: 1129
## simulation: 1130
## simulation: 1131
## simulation: 1132
## simulation: 1133
## simulation: 1134
## simulation: 1135
## simulation: 1136
## simulation: 1137
## simulation: 1138
## simulation: 1139
## simulation: 1140
## simulation: 1141
## simulation: 1142
## simulation: 1143
## simulation: 1144
## simulation: 1145
## simulation: 1146
## simulation: 1147
## simulation: 1148
## simulation: 1149
## simulation: 1150
## simulation: 1151
## simulation: 1152
## simulation: 1153
## simulation: 1154
## simulation: 1155
## simulation: 1156
## simulation: 1157
## simulation: 1158
## simulation: 1159
## simulation: 1160
## simulation: 1161
## simulation: 1162
## simulation: 1163
## simulation: 1164
## simulation: 1165
## simulation: 1166
## simulation: 1167
## simulation: 1168
## simulation: 1169
## simulation: 1170
## simulation: 1171
## simulation: 1172
## simulation: 1173
## simulation: 1174
## simulation: 1175
## simulation: 1176
## simulation: 1177
## simulation: 1178
## simulation: 1179
## simulation: 1180
## simulation: 1181
## simulation: 1182
## simulation: 1183
## simulation: 1184
## simulation: 1185
## simulation: 1186
## simulation: 1187
## simulation: 1188
## simulation: 1189
## simulation: 1190
## simulation: 1191
## simulation: 1192
## simulation: 1193
## simulation: 1194
## simulation: 1195
## simulation: 1196
## simulation: 1197
## simulation: 1198
## simulation: 1199
## simulation: 1200
Check sim data
head(sim_datasets)
## # A tibble: 6 x 12
## sim_num alpha0_mu alpha0_sigma alphaD_mu alphaD_sigma beta0_mu
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 3.65 0.165 -0.0314 0.141 -3.93
## 2 2 3.46 0.146 -0.447 0.122 2.01
## 3 3 4.32 0.417 -0.122 0.260 1.06
## 4 4 3.45 0.0352 0.161 0.000496 0.315
## 5 5 3.85 0.0252 -0.523 0.0131 0.0268
## 6 6 3.80 0.332 0.156 0.460 3.77
## # … with 6 more variables: beta0_sigma <dbl>, betaD_mu <dbl>,
## # betaD_sigma <dbl>, nsubj <dbl>, nobs_per_cond <dbl>, dataset <list>
head(sim_datasets$dataset[[1]])
## # A tibble: 6 x 7
## subj subj_alpha0 subj_alphaD subj_beta0 subj_betaD nobs_per_condit…
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 3.56 0.304 -4.20 -1.02 126
## 2 2 3.69 0.0279 -3.88 -0.173 126
## 3 3 3.83 -0.0175 -3.84 -0.809 126
## 4 4 3.75 -0.105 -4.32 -0.418 126
## 5 5 3.55 0.0680 -3.48 1.55 126
## 6 6 3.31 -0.118 -3.88 1.57 126
## # … with 1 more variable: subj_obs <list>
head(sim_datasets$dataset[[1]]$subj_obs[[1]])
## obs_radian obs_degree stimulation
## 1 0.04073781 2.334104 0
## 2 -0.40981437 -23.480634 0
## 3 2.23162565 127.862731 0
## 4 2.81399587 161.230087 0
## 5 1.07569022 61.632510 0
## 6 2.94621439 168.805650 0
Ideas:
Definitely plot the priors on the intercept and slope means, transformed back into sd
plot the prior predicitive densities for subjects
Some other plot ideas:
shaded histogram plot, combining data across stimulation conditions - check if data looks too peaked or too flat
shaded histogram plot, one for each stimulation condition - same check ^
plot of density lines, one from each sim dataset, one with both conditions and one with them split out - this would help identify weird distribution more that shaded histogram perhaps
sim_datasets %>%
select(alpha0_mu, alpha0_sigma, alphaD_mu, alphaD_sigma) %>%
ggpairs(progress = FALSE)
Good.
sim_datasets %>%
select(beta0_mu, beta0_sigma, betaD_mu, betaD_sigma) %>%
ggpairs(progress = FALSE)
Good.
params <- cross_df(list(pMem = seq(0, 1, 0.1), sd = seq(5, 18, 3)))
params %>%
mutate(sims = map2(params$pMem, sd2k_vec(pracma::deg2rad(params$sd)), ~simulateData_likelihood(.x, .y, 1000))) %>%
unnest(sims) %>%
ggplot(aes(x = obs_degree)) +
geom_histogram(binwidth = 1) +
facet_grid(rows = vars(pMem), cols = vars(sd), labeller = 'label_both', scales = 'free')
params <- cross_df(list(pMem = seq(0, 1, 0.1), sd = seq(20, 50, 5)))
params %>%
mutate(sims = map2(params$pMem, sd2k_vec(pracma::deg2rad(params$sd)), ~simulateData_likelihood(.x, .y, 1000))) %>%
unnest(sims) %>%
ggplot(aes(x = obs_degree)) +
geom_histogram(binwidth = 1) +
facet_grid(rows = vars(pMem), cols = vars(sd), labeller = 'label_both', scales = 'free')
params <- cross_df(list(pMem = seq(0, 1, 0.1), sd = seq(55, 105, 5)))
params %>%
mutate(sims = map2(params$pMem, sd2k_vec(pracma::deg2rad(params$sd)), ~simulateData_likelihood(.x, .y, 1000))) %>%
unnest(sims) %>%
ggplot(aes(x = obs_degree)) +
geom_histogram(binwidth = 1)+
facet_grid(rows = vars(pMem), cols = vars(sd), labeller = 'label_both')
(pre group mean, post group mean, post-pre mean change)
No accounting for subject variance here
##
## alpha0_mu
##
## linked prior distribution
alpha0_mu_hist <- sim_datasets %>%
ggplot(aes(x = alpha0_mu)) +
geom_histogram(binwidth = 0.05, fill = "#F16C66") +
theme(legend.position = "none",
axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
labs(x = "linked prior on group mean for circSD, pre condition",
title = "circSD group mean prior, alpha0_mu") +
scale_x_continuous(limits = c(1,7))
##
## alpha0_mu + alphaD_mu = alpha1_mu
##
## linked prior distribution
alpha1_mu_hist <- sim_datasets %>%
transmute(alpha1_mu = alpha0_mu + alphaD_mu) %>%
ggplot(aes(x = alpha1_mu)) +
geom_histogram(binwidth = 0.05, fill = "#685369") +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
labs(x = "linked prior on group mean for circSD, post condition",
title = "circSD group mean prior, alpha0_mu + alphaD_mu") +
scale_x_continuous(limits = c(1,7))
plot_grid(alpha0_mu_hist, alpha1_mu_hist, nrow = 1, align = "h")
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
##
## alpha0_mu
##
## delinked prior alpha0_mu
alpha0_mu_delinked_stats <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu)) %>%
summarise(less_than_5 = sum(delinked < 5)/n(), greater_than_120 = sum(delinked > 120)/n())
alpha0_mu_delinked_hist <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 2, fill = "#F16C66") +
geom_vline(xintercept = c(5, 120), linetype = "dashed") +
labs(x = "prior on group mean for circSD, pre condition",
title = "circSD group mean prior, exp(alpha0_mu)",
subtitle = glue("prob(circSD < 5) = {alpha0_mu_delinked_stats$less_than_5}\nprob(circSD > 120) = {alpha0_mu_delinked_stats$greater_than_120}")) +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12))
##
## alpha0_mu + alphaD_mu = alpha1_mu
##
## delinked prior alpha0_mu + alphaD_mu
alpha1_mu_delinked_hist <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu + alphaD_mu)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 2, fill = "#685369") +
labs(x = "prior on group mean for circSD, post condition",
title = "circSD group mean prior, exp(alpha0_mu + alphaD_mu)") +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12))
plot_grid(alpha0_mu_delinked_hist, alpha1_mu_delinked_hist, nrow = 1, align = "h")
Better containment of pre condition
##
## exp(alpha0_mu + alphaD_mu) - exp(alpha0_mu) plot
##
alphaD_mu_delinked_stats <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu + alphaD_mu) - exp(alpha0_mu)) %>%
summarise(less_than_n100 = sum(delinked < -100)/n(), greater_than_100 = sum(delinked > 100)/n())
alphaD_mu_delinked_hist <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu + alphaD_mu) - exp(alpha0_mu)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 2, fill = "#00AEB2") +
geom_vline(xintercept = c(-100, 100), linetype = "dashed") +
labs(x = "prior on group mean for delta-circSD, mean post minus mean pre",
title = "delta-circSD group mean prior, exp(alpha0_mu + alphaD_mu) - exp(alpha0_mu)",
subtitle = glue("prob(delta-circSD < -100) = {alphaD_mu_delinked_stats$less_than_n100}\nprob(delta-circSD > 100) = {alphaD_mu_delinked_stats$greater_than_100}")) +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
scale_x_continuous(breaks= pretty_breaks(10))
alphaD_mu_delinked_hist
mostly unchanged from m2
(based on simulated subjects from each dataset)
These plots indicate the effects of the priors on alpha0_mu, alpha0_sigma, alphaD_mu, alphaD_sigma.
(ignoring simulation sample sizes)
sim_datasets_unnest <- sim_datasets %>%
unnest(dataset) %>%
mutate(alpha0_mu_delinked = exp(alpha0_mu),
alpha1_mu_delinked = exp(alpha0_mu + alphaD_mu),
subj_pre_circSD = exp(subj_alpha0),
subj_post_circSD = exp(subj_alphaD + subj_alpha0),
subj_circSD_ES = subj_post_circSD - subj_pre_circSD)
subj_pre_circSD_plot <- sim_datasets_unnest %>%
mutate(subj_pre_circSD = if_else(subj_pre_circSD > 120, 120, subj_pre_circSD)) %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_pre_circSD)) +
geom_histogram(binwidth = 5, fill = "red", alpha = 0.3) +
labs(x = "subj_pre_circSD [prior predictive distribution]",
subtitle = glue("p(<5) = {sum(sim_datasets_unnest$subj_pre_circSD < 5)/length(sim_datasets_unnest$subj_pre_circSD)}\n p(>120) = {sum(sim_datasets_unnest$subj_pre_circSD > 120)/length(sim_datasets_unnest$subj_pre_circSD)}"))
subj_post_circSD_plot <- sim_datasets_unnest %>%
mutate(subj_post_circSD = if_else(subj_post_circSD > 120, 120, subj_post_circSD)) %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_post_circSD)) +
geom_histogram(binwidth = 5, fill = "red", alpha = 0.3) +
labs(x = "subj_post_circSD [prior predictive distribution]",
subtitle = glue("p(<5) = {sum(sim_datasets_unnest$subj_post_circSD < 5)/length(sim_datasets_unnest$subj_post_circSD)}\n p(>120) = {sum(sim_datasets_unnest$subj_post_circSD > 120)/length(sim_datasets_unnest$subj_post_circSD)}"))
subj_circSD_ES_plot <- sim_datasets_unnest %>%
mutate(subj_circSD_ES = if_else(subj_circSD_ES > 100, 100, subj_circSD_ES)) %>%
mutate(subj_circSD_ES = if_else(subj_circSD_ES < -100, -100, subj_circSD_ES)) %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_circSD_ES)) +
geom_histogram(binwidth = 5, fill = "red", alpha = 0.7) +
labs(x = "subj_circSD_ES [prior predictive distribution]",
subtitle = glue("p(< -100) = {sum(sim_datasets_unnest$subj_circSD_ES < -100)/length(sim_datasets_unnest$subj_circSD_ES)}\n p(>100) = {sum(sim_datasets_unnest$subj_circSD_ES > 100)/length(sim_datasets_unnest$subj_circSD_ES)}"))
plot_grid(subj_pre_circSD_plot, subj_post_circSD_plot, subj_circSD_ES_plot, ncol = 1)
Even better improvements in containment compared to m2.
Calculate min and max subject pre, post, and circSD effect size in each simulation.
sim_subj_extremes <- sim_datasets %>%
unnest(dataset) %>%
mutate(subj_alpha1 = subj_alpha0 + subj_alphaD,
subj_alpha0_delinked = exp(subj_alpha0),
subj_alpha1_delinked = exp(subj_alpha1),
subj_circSD_effect = subj_alpha1_delinked - subj_alpha0_delinked ) %>%
group_by(sim_num) %>%
summarize_at(vars(subj_alpha0_delinked, subj_alpha1_delinked, subj_circSD_effect), list(max = max, min = min))
simSD <- sim_datasets %>%
unnest(dataset) %>%
mutate(subj_pre_circSD = exp(subj_alpha0),
subj_post_circSD = exp(subj_alpha0 + subj_alphaD)) %>%
group_by(sim_num) %>%
summarise_at(vars(subj_pre_circSD, subj_post_circSD), list(mean = mean, sd = sd))
plot_grid(
ggplot(simSD, aes(subj_pre_circSD_sd)) +
geom_histogram(binwidth = 2) +
xlim(0, 200) +
labs(x = "pre condition: observed SD of distribution of circSD, per sim",
subtitle = glue::glue("p(>50) = {mean(simSD$subj_pre_circSD_sd > 50)}")),
ggplot(simSD, aes(subj_post_circSD_sd)) +
geom_histogram(binwidth = 2) +
xlim(0, 200) +
labs(x = "post condition: observed SD of distribution of circSD, per sim",
subtitle = glue::glue("p(>50) = {mean(simSD$subj_post_circSD_sd > 50)}")),
ncol = 1
)
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
plot_grid(
ggplot(simSD, aes(subj_pre_circSD_mean, subj_pre_circSD_sd)) +
geom_point() +
xlim(0, 200) + ylim(0, 100) +
labs(x = "pre condition: observed mean of distribution of circSD, per sim",
y = "pre condition: observed sd"),
ggplot(simSD, aes(subj_post_circSD_mean, subj_post_circSD_sd)) +
geom_point() +
xlim(0, 200) + ylim(0, 100) +
labs(x = "post condition: observed mean of distribution of circSD, per sim",
y = "post condition: observed sd"),
ncol = 1
)
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 35 rows containing missing values (geom_point).
#What is the max subject pre condition value in each dataset
pre_plot <- sim_subj_extremes %>%
mutate(subj_alpha0_delinked_max = if_else(subj_alpha0_delinked_max > 120, 120, subj_alpha0_delinked_max)) %>%
ggplot() +
geom_histogram(aes(x = subj_alpha0_delinked_max, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "pre condition: max subj circSD per sim [subj_alpha0_delinked_max]",
subtitle = glue("p(<5) = {sum(sim_subj_extremes$subj_alpha0_delinked_max < 5)/length(sim_subj_extremes$subj_alpha0_delinked_max)}\n p(>120) = {sum(sim_subj_extremes$subj_alpha0_delinked_max > 120)/length(sim_subj_extremes$subj_alpha0_delinked_max)}"))
#What is the max subject post condition value in each dataset
post_plot <- sim_subj_extremes %>%
mutate(subj_alpha1_delinked_max = if_else(subj_alpha1_delinked_max > 120, 120, subj_alpha1_delinked_max)) %>%
ggplot() +
geom_histogram(aes(x = subj_alpha1_delinked_max, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "post condition: max subj circSD per sim [subj_alpha1_delinked_max]",
subtitle = glue("p(<5) = {sum(sim_subj_extremes$subj_alpha1_delinked_max < 5)/length(sim_subj_extremes$subj_alpha1_delinked_max)}\n p(>120) = {sum(sim_subj_extremes$subj_alpha1_delinked_max > 120)/length(sim_subj_extremes$subj_alpha1_delinked_max)}"))
plot_grid(pre_plot, post_plot, ncol =1, align = 'v')
better containment than m2. m4 should do better by narrowing hierarchical mean priors.
pre_plot <- sim_subj_extremes %>%
mutate(subj_alpha0_delinked_min = if_else(subj_alpha0_delinked_min > 120, 120, subj_alpha0_delinked_min)) %>%
ggplot() +
geom_histogram(aes(x = subj_alpha0_delinked_min, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "pre condition: min subj circSD per sim [subj_alpha0_delinked_min]",
subtitle = glue("p(<5) = {sum(sim_subj_extremes$subj_alpha0_delinked_min < 5)/length(sim_subj_extremes$subj_alpha0_delinked_min)}\n p(>120) = {sum(sim_subj_extremes$subj_alpha0_delinked_min > 120)/length(sim_subj_extremes$subj_alpha0_delinked_min)}"))
#What is the most extreme subject post condition value in each dataset
post_plot <- sim_subj_extremes %>%
mutate(subj_alpha1_delinked_min = if_else(subj_alpha1_delinked_min > 120, 120, subj_alpha1_delinked_min)) %>%
ggplot() +
geom_histogram(aes(x = subj_alpha1_delinked_min, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "post condition: min subj circSD per sim [subj_alpha1_delinked_min]",
subtitle = glue("p(<5) = {sum(sim_subj_extremes$subj_alpha1_delinked_min < 5)/length(sim_subj_extremes$subj_alpha1_delinked_min)}\n p(>120) = {sum(sim_subj_extremes$subj_alpha1_delinked_min > 120)/length(sim_subj_extremes$subj_alpha1_delinked_min)}"))
plot_grid(pre_plot, post_plot, ncol =1, align = 'v')
Things look OK from this point of view.
#What is the most extreme change from pre to post in a subject
min_plot <- sim_subj_extremes %>%
mutate(subj_circSD_effect_min = if_else(subj_circSD_effect_min > 100 , 100, subj_circSD_effect_min)) %>%
mutate(subj_circSD_effect_min = if_else(subj_circSD_effect_min < -100 , -100, subj_circSD_effect_min)) %>%
ggplot() +
geom_histogram(aes(x = subj_circSD_effect_min, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "min subj circSD effect per sim \n[min(subj_alpha1_delinked - subj_alpha0_delinked)]",
subtitle = glue("p(<-100) = {sum(sim_subj_extremes$subj_circSD_effect_min < -100)/length(sim_subj_extremes$subj_circSD_effect_min)}\n p(>100) = {sum(sim_subj_extremes$subj_circSD_effect_min > 100)/length(sim_subj_extremes$subj_circSD_effect_min)}"))
#What is the most extreme subject post condition value in each dataset
max_plot <- sim_subj_extremes %>%
mutate(subj_circSD_effect_max = if_else(subj_circSD_effect_max > 100 , 100, subj_circSD_effect_max)) %>%
mutate(subj_circSD_effect_max = if_else(subj_circSD_effect_max < -100 , -100, subj_circSD_effect_max)) %>%
ggplot() +
geom_histogram(aes(x = subj_circSD_effect_max, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "max subj circSD effect per sim \n[max(subj_alpha1_delinked - subj_alpha0_delinked)]",
subtitle = glue("p(<-100) = {sum(sim_subj_extremes$subj_circSD_effect_max < -100)/length(sim_subj_extremes$subj_circSD_effect_max)}\n p(>100) = {sum(sim_subj_extremes$subj_circSD_effect_max > 100)/length(sim_subj_extremes$subj_circSD_effect_max)}"))
plot_grid(min_plot, max_plot, ncol =1, align = 'v')
maximum effects are more contained, very slightly. not a huge improvement there.
Both the max pre/post subject distributions of circSD and the max/min subject effect size distributions on circSD contained too many extreme values.
Can I see which prior distributions I need to change in order to correct these by plotting the points from those distributions along with their respective prior distribution samples?
(Is this a useful technique in general?)
Plots:
for (max pre circSD, max post circSD, max circSD ES, min circSD ES) * vs alpha0_mu * vs alpha0_sigma * vs alphaD_mu * vs alphaD_sigma
sim_datasets_join_extremes <- sim_datasets %>%
mutate(alpha0_mu_delinked = exp(alpha0_mu),
alpha1_mu_delinked = exp(alpha0_mu + alphaD_mu)) %>%
left_join(sim_subj_extremes, by = "sim_num") %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup()
#parameters describing uncertanity over group circ SDs in the pre condition (alpha0_mu, alpha0_sigma)
#vs
#implied pre condition min and max circSDs
preMax_vs_alpha0_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha0_delinked_max, y = alpha0_mu_delinked)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
preMin_vs_alpha0_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha0_delinked_min, y = alpha0_mu_delinked)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
preMax_vs_alpha0_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha0_delinked_max, y = alpha0_sigma)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
preMin_vs_alpha0_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha0_delinked_min, y = alpha0_sigma)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
plot_grid(preMin_vs_alpha0_mu, preMax_vs_alpha0_mu, preMin_vs_alpha0_sigma, preMax_vs_alpha0_sigma, nrow = 2, align = "hv")
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 25 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 25 rows containing missing values (geom_point).
#parameters describing uncertanity over change from pre to post condition (alphaD_mu, alphaD_sigma)
#vs
#implied post condition min and max circSDs
postMin_vs_alphaD_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha1_delinked_min, y = alphaD_mu)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
postMax_vs_alphaD_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha1_delinked_max, y = alphaD_mu)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
postMin_vs_alphaD_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha1_delinked_min, y = alphaD_sigma)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
postMax_vs_alphaD_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha1_delinked_max, y = alphaD_sigma)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
plot_grid(postMin_vs_alphaD_mu, postMax_vs_alphaD_mu, postMin_vs_alphaD_sigma, postMax_vs_alphaD_sigma, nrow = 2, align = "hv")
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 119 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 119 rows containing missing values (geom_point).
#parameters describing uncertanity over group circ SDs in the pre condition (alpha0_mu, alpha0_sigma)
#vs
#implied effect size min and max
ESmin_vs_alpha0_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_min, y = alpha0_mu_delinked)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmax_vs_alpha0_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_max, y = alpha0_mu_delinked)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmin_vs_alpha0_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_min, y = alpha0_sigma)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmax_vs_alpha0_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_max, y = alpha0_sigma)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
plot_grid(ESmin_vs_alpha0_mu, ESmax_vs_alpha0_mu, ESmin_vs_alpha0_sigma, ESmax_vs_alpha0_sigma, nrow = 2, align = "hv")
## Warning: Removed 19 rows containing missing values (geom_point).
## Warning: Removed 139 rows containing missing values (geom_point).
## Warning: Removed 19 rows containing missing values (geom_point).
## Warning: Removed 139 rows containing missing values (geom_point).
#parameters describing uncertanity over change from pre to post condition(alphaD_mu, alphaD_sigma)
#vs
#implied effect size min and max
ESmin_vs_alphaD_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_min, y = alphaD_mu)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmax_vs_alphaD_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_max, y = alphaD_mu)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmin_vs_alphaD_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_min, y = alphaD_sigma)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmax_vs_alphaD_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_max, y = alphaD_sigma)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
plot_grid(ESmin_vs_alphaD_mu, ESmax_vs_alphaD_mu, ESmin_vs_alphaD_sigma, ESmax_vs_alphaD_sigma, nrow = 2, align = "hv")
## Warning: Removed 19 rows containing missing values (geom_point).
## Warning: Removed 139 rows containing missing values (geom_point).
## Warning: Removed 19 rows containing missing values (geom_point).
## Warning: Removed 139 rows containing missing values (geom_point).
Takeaways:
Unsurprisingly, effect size and pre,post extreme values would likely be better contained with narrower priors.
First ideas:
*Only narrow priors on mu parameters, alpha0_mu and alphaD_mu. This should help but the next step will likely be ncessary?
*Also narrow priors on the SD parameters, alpha0_sigma and alphaD_sigma.
(pre group mean, post group mean, post-pre mean change)
## linked prior distribution
plot_grid(
sim_datasets %>%
ggplot(aes(x = beta0_mu)) +
geom_histogram(binwidth = 0.05, fill = "#F16C66") +
theme(legend.position = "none",
axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
labs(x = "linked prior on group mean for pMem, pre condition",
title = "pMem group mean prior, beta0_mu") +
scale_x_continuous(limits = c(-5,5))
,
sim_datasets %>%
transmute(beta1_mu = beta0_mu + betaD_mu) %>%
ggplot(aes(x = beta1_mu)) +
geom_histogram(binwidth = 0.05, fill = "#685369") +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
labs(x = "linked prior on group mean for pMem, post condition",
title = "pMem group mean prior, beta0_mu + betaD_mu") +
scale_x_continuous(limits = c(-5,5))
,
nrow = 1,
align = "h"
)
## Warning: Removed 1 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 4 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
##
## inv_logit(beta0_mu)
##
##
## inv_logit(beta0_mu + betaD_mu) = inv_logit(beta1_mu)
##
## delinked prior alpha0_mu
beta0_mu_delinked_stats <- sim_datasets %>%
transmute(delinked = exp(beta0_mu)/(exp(beta0_mu) + 1)) %>%
summarise(less_than_5 = sum(delinked < 0.05)/n(), greater_than_95 = sum(delinked > 0.95)/n())
plot_grid(
sim_datasets %>%
transmute(delinked = exp(beta0_mu)/(exp(beta0_mu) + 1)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 0.03, fill = "#F16C66") +
geom_vline(xintercept = c(0.05, 0.95), linetype = "dashed") +
labs(x = "prior on group mean for pMem, pre condition",
title = "pMem group mean prior, inv_logit(beta0_mu)",
subtitle = glue("prob(pMem < 0.05) = {beta0_mu_delinked_stats$less_than_5}\nprob(pMem > 0.95) = {beta0_mu_delinked_stats$greater_than_95}")) +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12))
,
sim_datasets %>%
transmute(delinked = exp(beta0_mu + betaD_mu)/(exp(beta0_mu + betaD_mu) + 1)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 0.03, fill = "#685369") +
labs(x = "prior on group mean for pMem, post condition",
title = "pMem group mean prior, inv_logit(alpha0_mu + alphaD_mu)") +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12))
,
nrow = 1,
align = "h"
)
##
## inv_logit(beta0_mu + betaD_mu) - inv_logit(beta0_mu) plot
##
betaD_mu_delinked_stats <- sim_datasets %>%
transmute(delinked = exp(beta0_mu + betaD_mu)/(exp(beta0_mu + betaD_mu) + 1) - exp(beta0_mu)/(exp(beta0_mu) + 1)) %>%
summarise(less_than_n80 = sum(delinked < -0.8)/n(), greater_than_80 = sum(delinked > 0.8)/n())
sim_datasets %>%
transmute(delinked = exp(beta0_mu + betaD_mu)/(exp(beta0_mu + betaD_mu) + 1) - exp(beta0_mu)/(exp(beta0_mu) + 1)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 0.02, fill = "#00AEB2") +
geom_vline(xintercept = c(-0.8, 0.8), linetype = "dashed") +
labs(x = "prior on group mean for delta-pMem, mean post minus mean pre",
title = "delta-pMem group mean prior, inv_logit(beta0_mu + betaD_mu) - inv_logit(beta0_mu)",
subtitle = glue("prob(delta-pMem < -0.8) = {betaD_mu_delinked_stats$less_than_n80}\nprob(delta-pMem > 0.8) = {betaD_mu_delinked_stats$greater_than_80}")) +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
scale_x_continuous(breaks= pretty_breaks(10))
(based on simulated subjects from each dataset)
(ignoring simulation sample sizes)
sim_datasets_unnest <- sim_datasets %>%
unnest(dataset) %>%
mutate(beta0_mu_delinked = exp(beta0_mu)/(exp(beta0_mu) + 1),
beta1_mu_delinked = exp(beta0_mu + betaD_mu)/(exp(beta0_mu + betaD_mu) + 1),
subj_pre_pMem = exp(subj_beta0)/(exp(subj_beta0) + 1),
subj_post_pMem = exp(subj_betaD + subj_beta0)/(exp(subj_betaD + subj_beta0) + 1),
subj_pMem_ES = subj_post_pMem - subj_pre_pMem)
pMem_ES_quantiles <- quantile(sim_datasets_unnest$subj_pMem_ES, c(0.05, 0.95))
plot_grid(
sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_pre_pMem)) +
geom_histogram(binwidth = 0.01, fill = "red", alpha = 0.3) +
labs(x = "subj_pre_pMem [prior predictive distribution]",
subtitle = glue("p(< 0.05) = {sum(sim_datasets_unnest$subj_pre_pMem < 0.05)/length(sim_datasets_unnest$subj_pre_pMem)}\n p(> 0.95) = {sum(sim_datasets_unnest$subj_pre_pMem > 0.95)/length(sim_datasets_unnest$subj_pre_pMem)}"))
,
sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_post_pMem)) +
geom_histogram(binwidth = 0.01, fill = "red", alpha = 0.3) +
labs(x = "subj_post_pMem [prior predictive distribution]",
subtitle = glue("p(< 0.05) = {sum(sim_datasets_unnest$subj_post_pMem < 0.05)/length(sim_datasets_unnest$subj_post_pMem)}\n p(> 0.95) = {sum(sim_datasets_unnest$subj_post_pMem > 0.95)/length(sim_datasets_unnest$subj_post_pMem)}"))
,
sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_pMem_ES)) +
geom_histogram(binwidth = 0.01, fill = "red", alpha = 0.7) +
geom_vline(xintercept = pMem_ES_quantiles, linetype = "dashed") +
labs(x = "subj_pMem_ES [prior predictive distribution], (5%, 95%) quantile lines",
subtitle = glue("p(< -0.8) = {sum(sim_datasets_unnest$subj_pMem_ES < -0.8)/length(sim_datasets_unnest$subj_pMem_ES)}\n p(> 0.8) = {sum(sim_datasets_unnest$subj_pMem_ES > 0.80)/length(sim_datasets_unnest$subj_pMem_ES)}"))
,
ncol = 1
)
Calculate min and max subject pre, post, and pMem effect size in each simulation.
sim_subj_extremes <- sim_datasets %>%
unnest(dataset) %>%
mutate(subj_beta1 = subj_beta0 + subj_betaD,
subj_beta0_delinked = exp(subj_beta0)/(exp(subj_beta0) + 1),
subj_beta1_delinked = exp(subj_beta0 + subj_betaD)/(exp(subj_beta0 + subj_betaD) + 1),
subj_pMem_effect = subj_beta1_delinked - subj_beta0_delinked ) %>%
group_by(sim_num) %>%
summarize_at(vars(subj_beta0_delinked, subj_beta1_delinked, subj_pMem_effect), list(max = max, min = min))
plot_grid(
#What is the max subject pre condition value in each dataset
sim_subj_extremes %>%
ggplot() +
geom_histogram(aes(x = subj_beta0_delinked_max, fill = "red", alpha = 0.3), binwidth = 0.02) +
theme(legend.position = "none") +
labs(x = "pre condition: max subj pMem per sim [subj_beta0_delinked_max]",
subtitle = glue("p(< 0.05) = {sum(sim_subj_extremes$subj_beta0_delinked_max < 0.05)/length(sim_subj_extremes$subj_beta0_delinked_max)}\n p(> 0.95) = {sum(sim_subj_extremes$subj_beta0_delinked_max > 0.95)/length(sim_subj_extremes$subj_beta0_delinked_max)}"))
,
#What is the max subject post condition value in each dataset
sim_subj_extremes %>%
ggplot() +
geom_histogram(aes(x = subj_beta1_delinked_max, fill = "red", alpha = 0.3), binwidth = 0.02) +
theme(legend.position = "none") +
labs(x = "post condition: max subj pMem per sim [subj_beta1_delinked_max]",
subtitle = glue("p(< 0.05) = {sum(sim_subj_extremes$subj_beta1_delinked_max < 0.05)/length(sim_subj_extremes$subj_beta1_delinked_max)}\n p(>0.95) = {sum(sim_subj_extremes$subj_beta1_delinked_max > 0.95)/length(sim_subj_extremes$subj_beta1_delinked_max)}"))
,
ncol =1,
align = 'v'
)
Perhaps too skewed
plot_grid(
#What is the min subject pre condition value in each dataset
sim_subj_extremes %>%
ggplot() +
geom_histogram(aes(x = subj_beta0_delinked_min, fill = "red", alpha = 0.3), binwidth = 0.02) +
theme(legend.position = "none") +
labs(x = "pre condition: min subj pMem per sim [subj_beta0_delinked_min]",
subtitle = glue("p(< 0.05) = {sum(sim_subj_extremes$subj_beta0_delinked_min < 0.05)/length(sim_subj_extremes$subj_beta0_delinked_min)}\n p(> 0.95) = {sum(sim_subj_extremes$subj_beta0_delinked_min > 0.95)/length(sim_subj_extremes$subj_beta0_delinked_min)}"))
,
#What is the min subject post condition value in each dataset
sim_subj_extremes %>%
ggplot() +
geom_histogram(aes(x = subj_beta1_delinked_min, fill = "red", alpha = 0.3), binwidth = 0.02) +
theme(legend.position = "none") +
labs(x = "post condition: min subj pMem per sim [subj_beta1_delinked_min]",
subtitle = glue("p(< 0.05) = {sum(sim_subj_extremes$subj_beta1_delinked_min < 0.05)/length(sim_subj_extremes$subj_beta1_delinked_min)}\n p(>0.95) = {sum(sim_subj_extremes$subj_beta1_delinked_min > 0.95)/length(sim_subj_extremes$subj_beta1_delinked_min)}"))
,
ncol =1,
align = 'v'
)
#What is the most extreme change from pre to post in a subject
min_plot <- sim_subj_extremes %>%
mutate(subj_pMem_effect_min = if_else(subj_pMem_effect_min > 1 , 1, subj_pMem_effect_min)) %>%
mutate(subj_pMem_effect_min = if_else(subj_pMem_effect_min < -1 , -1, subj_pMem_effect_min)) %>%
ggplot() +
geom_histogram(aes(x = subj_pMem_effect_min, fill = "red", alpha = 0.3), binwidth = 0.03) +
theme(legend.position = "none") +
labs(x = "min subj pMem effect per sim \n[min(subj_beta1_delinked - subj_beta0_delinked)]",
subtitle = glue("p(<-0.8) = {sum(sim_subj_extremes$subj_pMem_effect_min < -0.8)/length(sim_subj_extremes$subj_pMem_effect_min)}\n p(>0.8) = {sum(sim_subj_extremes$subj_pMem_effect_min > 0.8)/length(sim_subj_extremes$subj_pMem_effect_min)}")) +
xlim(c(-1, 1))
#What is the most extreme subject post condition value in each dataset
max_plot <- sim_subj_extremes %>%
mutate(subj_pMem_effect_max = if_else(subj_pMem_effect_max > 1 , 1, subj_pMem_effect_max)) %>%
mutate(subj_pMem_effect_max = if_else(subj_pMem_effect_max < -1 , -1, subj_pMem_effect_max)) %>%
ggplot() +
geom_histogram(aes(x = subj_pMem_effect_max, fill = "red", alpha = 0.3), binwidth = 0.03) +
theme(legend.position = "none") +
labs(x = "max subj pMem effect per sim \n[max(subj_beta1_delinked - subj_beta0_delinked)]",
subtitle = glue("p(<-0.8) = {sum(sim_subj_extremes$subj_pMem_effect_max < -0.8)/length(sim_subj_extremes$subj_pMem_effect_max)}\n p(>0.8) = {sum(sim_subj_extremes$subj_pMem_effect_max > 0.8)/length(sim_subj_extremes$subj_pMem_effect_max)}")) +
xlim(c(-1, 1))
plot_grid(min_plot, max_plot, ncol =1, align = 'v')
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
This looks good, could even be wider.
joint_pre_plot <- sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(y = exp(subj_alpha0), x = exp(subj_beta0)/(exp(subj_beta0) + 1))) +
geom_point() +
labs(x = 'pMem, pre', y = 'circSD, pre')
plot_grid(
joint_pre_plot,
ncol = 1,
align = 'v'
)
joint_ES_plot <- sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
mutate(subj_pMem_ES = exp(subj_beta0 + subj_betaD)/(exp(subj_beta0 + subj_betaD) + 1) - exp(subj_beta0)/(exp(subj_beta0) + 1),
subj_circSD_ES = exp(subj_alpha0 + subj_alphaD) - exp(subj_alpha0)) %>%
ggplot(aes(x = subj_pMem_ES, y = subj_circSD_ES)) +
geom_point() +
geom_vline(xintercept = 0, linetype = 'dashed') +
geom_hline(yintercept = 0, linetype = 'dashed')
plot_grid(
joint_ES_plot,
joint_ES_plot + ylim(c(-200, 200)) + labs(subtitle = "zoom in"),
ncol = 1,
align = 'v'
)
## Warning: Removed 10 rows containing missing values (geom_point).
(based on simulated obs from each dataset, ignoring subject labels)
sim_subj_obs_hist_count <- function(dataset, condition = 0){
dataset_obs <- dataset %>%
unnest(subj_obs) %>%
filter(stimulation == condition) %>%
select(obs_degree)
breaks <- seq(-180, 180, 5)
bincount <- hist(dataset_obs$obs_degree, breaks = breaks, plot = FALSE)$counts
bincount_names <- glue("c{breaks[-1]}")
names(bincount) <- bincount_names
bincount_df <- data.frame(as.list(bincount))
return(bincount_df)
}
make_quantmat <- function(sim_datasets, condition = 0){
bincounts <- sim_datasets %>%
select(dataset) %>%
mutate(subj_hist_counts = map(dataset, sim_subj_obs_hist_count, condition)) %>%
select(-dataset) %>%
unnest(subj_hist_counts) %>%
as_tibble()
xvals <- seq(-177.5, 177.5, 5)
probs <- seq(0.1,0.9,0.1)
quantmat <- as.data.frame(matrix(NA, nrow=ncol(bincounts), ncol=length(probs)))
names(quantmat) <- paste0("p",probs)
quantmat <- cbind(quantmat, xvals)
for (iQuant in 1:length(probs)){
quantmat[,paste0("p",probs[iQuant])] <- as.numeric(summarise_all(bincounts, ~quantile(., probs[iQuant])))
}
return(quantmat)
}
# calculate quantile mats from each condition
quantmat_cond0 <- make_quantmat(sim_datasets, 0)
quantmat_cond1 <- make_quantmat(sim_datasets, 1)
c_light <- "#DCBCBC"
c_light_highlight <- "#C79999"
c_mid <- "#B97C7C"
c_mid_highlight <- "#A25050"
c_dark <- "#8F2727"
c_dark_highlight <- "#7C0000"
plot_grid(
ggplot(quantmat_cond0, aes(x = xvals)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), color = c_dark, size = 1) +
scale_x_continuous(breaks=pretty_breaks(10)) +
labs(x = "error (degrees)", y = "count +/- quantile", subtitle = "without stimulation")
,
ggplot(quantmat_cond1, aes(x = xvals)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), color = c_dark, size = 1) +
scale_x_continuous(breaks=pretty_breaks(10)) +
labs(x = "error (degrees)", y = "count +/- quantile", subtitle = "with stimulation")
,
sim_datasets %>%
sample_n(1) %>%
unnest(dataset) %>%
unnest(subj_obs) %>%
filter(stimulation == sample(c(0,1), 1)) %T>%
print() %>%
ggplot(aes(x = obs_degree)) +
geom_density() +
labs(subtitle = "correctness check: random simulation, random condition") +
scale_x_continuous(breaks=pretty_breaks(10))
,
ncol = 1,
align = "v"
)
## # A tibble: 1,890 x 20
## sim_num alpha0_mu alpha0_sigma alphaD_mu alphaD_sigma beta0_mu
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1019 4.07 0.256 0.333 0.383 -0.324
## 2 1019 4.07 0.256 0.333 0.383 -0.324
## 3 1019 4.07 0.256 0.333 0.383 -0.324
## 4 1019 4.07 0.256 0.333 0.383 -0.324
## 5 1019 4.07 0.256 0.333 0.383 -0.324
## 6 1019 4.07 0.256 0.333 0.383 -0.324
## 7 1019 4.07 0.256 0.333 0.383 -0.324
## 8 1019 4.07 0.256 0.333 0.383 -0.324
## 9 1019 4.07 0.256 0.333 0.383 -0.324
## 10 1019 4.07 0.256 0.333 0.383 -0.324
## # … with 1,880 more rows, and 14 more variables: beta0_sigma <dbl>,
## # betaD_mu <dbl>, betaD_sigma <dbl>, nsubj <dbl>, nobs_per_cond <dbl>,
## # subj <int>, subj_alpha0 <dbl>, subj_alphaD <dbl>, subj_beta0 <dbl>,
## # subj_betaD <dbl>, nobs_per_condition <dbl>, obs_radian <dbl>,
## # obs_degree <dbl>, stimulation <dbl>
looks good.